General Information

sapply(dat, summary) %>% t() %>% kable(caption = "General summary of Boston dataset variables") %>% 
    kable_styling(position = "left", bootstrap_options = c("striped", "hover")) %>% 
    column_spec(1, bold = T)

(a) Simple linear regression for each variable

# run regression
nm <- names(dat)[!grepl("crim", names(dat))]
fit.simple <- lapply(nm, function(x) {
    lm(as.formula(paste0("crim ~ ", x)), data = dat)
})

# print lm summary, graphs
invisible(lapply(fit.simple, mySummary))

Model: crim = 4.454 × (Intercept) + -0.074 × zn

Model: crim = -2.064 × (Intercept) + 0.51 × indus

Model: crim = 3.744 × (Intercept) + -1.893 × chas

Model: crim = -13.72 × (Intercept) + 31.249 × nox

Model: crim = 20.482 × (Intercept) + -2.684 × rm

Model: crim = -3.778 × (Intercept) + 0.108 × age

Model: crim = 9.499 × (Intercept) + -1.551 × dis

Model: crim = -2.287 × (Intercept) + 0.618 × rad

Model: crim = -8.528 × (Intercept) + 0.03 × tax

Model: crim = -17.647 × (Intercept) + 1.152 × ptratio

Model: crim = 16.554 × (Intercept) + -0.036 × black

Model: crim = -3.331 × (Intercept) + 0.549 × lstat

Model: crim = 11.797 × (Intercept) + -0.363 × medv

Statistically significant association between the predictor and the response variable has been observed for all the variables except chas.

(b) Multiple regression with all the variables

# run regression
fit.all <- lm(crim ~ ., data = dat)
mySummary(fit.all)

Model: crim = 17.033 × (Intercept) + 0.045 × zn + -0.064 × indus + -0.749 × chas + -10.314 × nox + 0.43 × rm + 0.001 × age + -0.987 × dis + 0.588 × rad + -0.004 × tax + -0.271 × ptratio + -0.008 × black + 0.126 × lstat + -0.199 × medv

Significant variables are: zn, dis, rad, black, medv.

(c) (a) and (b) results comparison plot

x <- sapply(fit.simple, function(x) coefficients(x)[2])
y <- coefficients(fit.all)[-1]
df <- data.table(Variable = names(x), Simple = x, Multiple = y, `Absolute difference` = abs(x - 
    y))

fig <- plot_ly(df, x = ~Simple, y = ~Multiple, type = "scatter", mode = "markers", 
    hoverinfo = "text", text = ~paste("</br> Variable: ", Variable, "</br> Simple regression: ", 
        round(Simple, 3), "</br> Multiple regression: ", round(Multiple, 3)))

fig <- fig %>% layout(title = "(a) and (b) results comparison plot", xaxis = list(title = "Simple regression for each variable"), 
    yaxis = list(title = "Multiple regression"))
fig

df[order(-`Absolute difference`)] %>% xtable() %>% kable() %>% kable_styling(position = "left", 
    bootstrap_options = "hover")

(d) Cubic model

# run regression
df.num <- Filter(function(x) {
    class(x) == "numeric"
}, dat)
nm.num <- names(df.num)[!grepl("crim", names(df.num))]
fit.cubic <- lapply(nm.num, function(x) {
    lm(as.formula(paste0("crim ~ ", paste0("poly(", x, ",3)"))), data = df.num)
})

# print lm summary, graphs
invisible(lapply(fit.cubic, mySummary))

Model: crim = 3.614 × (Intercept) + -38.75 × poly(zn, 3)1 + 23.94 × poly(zn, 3)2 + -10.072 × poly(zn, 3)3

Model: crim = 3.614 × (Intercept) + 78.591 × poly(indus, 3)1 + -24.395 × poly(indus, 3)2 + -54.13 × poly(indus, 3)3

Model: crim = 3.614 × (Intercept) + 81.372 × poly(nox, 3)1 + -28.829 × poly(nox, 3)2 + -60.362 × poly(nox, 3)3

Model: crim = 3.614 × (Intercept) + -42.379 × poly(rm, 3)1 + 26.577 × poly(rm, 3)2 + -5.51 × poly(rm, 3)3

Model: crim = 3.614 × (Intercept) + 68.182 × poly(age, 3)1 + 37.484 × poly(age, 3)2 + 21.353 × poly(age, 3)3

Model: crim = 3.614 × (Intercept) + -73.389 × poly(dis, 3)1 + 56.373 × poly(dis, 3)2 + -42.622 × poly(dis, 3)3

Model: crim = 3.614 × (Intercept) + 112.646 × poly(tax, 3)1 + 32.087 × poly(tax, 3)2 + -7.997 × poly(tax, 3)3

Model: crim = 3.614 × (Intercept) + 56.045 × poly(ptratio, 3)1 + 24.775 × poly(ptratio, 3)2 + -22.28 × poly(ptratio, 3)3

Model: crim = 3.614 × (Intercept) + -74.431 × poly(black, 3)1 + 5.926 × poly(black, 3)2 + -4.835 × poly(black, 3)3

Model: crim = 3.614 × (Intercept) + 88.07 × poly(lstat, 3)1 + 15.888 × poly(lstat, 3)2 + -11.574 × poly(lstat, 3)3

Model: crim = 3.614 × (Intercept) + -75.058 × poly(medv, 3)1 + 88.086 × poly(medv, 3)2 + -48.033 × poly(medv, 3)3

(a) Create train and test sets

set.seed(1995)
train_ind <- sample(nrow(OJ), 800)
train <- as.data.table(OJ[train_ind, ])
test <- as.data.table(OJ[-train_ind, ])

(b) Fit a tree

fit.tree = tree(Purchase ~ ., data = train)
summary(fit.tree)

Classification tree:
tree(formula = Purchase ~ ., data = train)
Variables actually used in tree construction:
[1] "LoyalCH"     "SpecialCH"   "SalePriceMM" "PriceDiff"  
Number of terminal nodes:  7 
Residual mean deviance:  0.7772 = 616.3 / 793 
Misclassification error rate: 0.1625 = 130 / 800 

4 variables were actually used in tree construction: LoyalCH, SpecialCH, SalePriceMM, PriceDiff.
Training error is 0.1462.
Tree has 7 terminal nodes.

(c) Tree info

fit.tree

node), split, n, deviance, yval, (yprob)
      * denotes terminal node

 1) root 800 1066.000 CH ( 0.61500 0.38500 )  
   2) LoyalCH < 0.450956 277  294.600 MM ( 0.22383 0.77617 )  
     4) LoyalCH < 0.0356415 53    9.922 MM ( 0.01887 0.98113 ) *
     5) LoyalCH > 0.0356415 224  262.300 MM ( 0.27232 0.72768 )  
      10) SpecialCH < 0.5 202  219.200 MM ( 0.23267 0.76733 )  
        20) SalePriceMM < 2.04 103   78.100 MM ( 0.12621 0.87379 ) *
        21) SalePriceMM > 2.04 99  127.400 MM ( 0.34343 0.65657 ) *
      11) SpecialCH > 0.5 22   28.840 CH ( 0.63636 0.36364 ) *
   3) LoyalCH > 0.450956 523  489.600 CH ( 0.82218 0.17782 )  
     6) LoyalCH < 0.705699 227  290.800 CH ( 0.66079 0.33921 )  
      12) PriceDiff < 0.085 81  107.800 MM ( 0.38272 0.61728 ) *
      13) PriceDiff > 0.085 146  139.800 CH ( 0.81507 0.18493 ) *
     7) LoyalCH > 0.705699 296  124.500 CH ( 0.94595 0.05405 ) *

Let’s pick the terminal node 11) SpecialCH > 0.5 22 29.77 CH ( 0.59091 0.40909 ) *.
This node is a child of 5) LoyalCH > 0.051325 232 269.40 MM ( 0.26724 0.73276 ) and the root node 2) LoyalCH < 0.48285 298 307.40 MM ( 0.21141 0.78859 ).
Split of this node is SpecialCH > 0.5, there are 22 observations (59% CH and 41% MM) in this node.
The final prediction for this branch is CH.

(d) Tree plot

plot(fit.tree, uniform = TRUE, main = "Classification Tree for Orange Juice Data Purchase")
text(fit.tree, cex = 0.9)

Only 4 variables out of 17 were used.
The most important variable in the tree is LoyalCH.
If LoyalCH < 0.051, the tree always predicts MM (Minute Maid) or if LoyalCH > 0.765, the tree always predicts CH (Citrus Hill).
In other ways (0.051 < LoyalCH < 0.765) decision also depends on SpecialCH or PriceDiff.

(e) Predictions

pred <- predict(fit.tree, test, type = "class")
conf.m <- caret::confusionMatrix(test$Purchase, pred)
test.acc <- round(as.numeric(conf.m$overall[1]), 3)

gg <- conf.m$table %>% data.frame() %>% mutate(Prediction = factor(Prediction, levels = c("MM", 
    "CH"))) %>% group_by(Reference) %>% mutate(total = sum(Freq)) %>% ungroup() %>% 
    ggplot(aes(Reference, Prediction, fill = Freq)) + geom_tile() + geom_text(aes(label = Freq), 
    size = 8) + scale_fill_gradient(low = "white", high = "#1a7ab1") + scale_x_discrete(position = "top") + 
    geom_tile(color = "black", fill = "black", alpha = 0)

gg

Accuracy on test set is: 0.796, then test error rate is: 0.204.

(f)

set.seed(1)
my.cv.tree <- cv.tree(fit.tree)

(g) and (h)

dd <- data.frame(size = my.cv.tree$size, dev = my.cv.tree$dev)

fig <- plot_ly(dd, x = ~size, y = ~dev, type = "scatter", mode = "lines", hoverinfo = "text", 
    text = ~paste("</br> Tree Size: ", size, "</br> Deviance: ", round(dev, 3)))

fig <- fig %>% layout(title = "Tree size vs deviance", xaxis = list(title = "Tree Size"), 
    yaxis = list(title = "Deviance"))
fig

Minimum deviance on cross validated data reached when tree size is 6.

(i)

fit.pruned <- prune.tree(fit.tree, best = 6)
plot(fit.pruned, uniform = TRUE, main = "Classification Tree for Orange Juice Data Purchase")
text(fit.pruned, cex = 0.9)

(j)

summary(fit.pruned)

Classification tree:
snip.tree(tree = fit.tree, nodes = 10L)
Variables actually used in tree construction:
[1] "LoyalCH"   "SpecialCH" "PriceDiff"
Number of terminal nodes:  6 
Residual mean deviance:  0.7935 = 630 / 794 
Misclassification error rate: 0.1625 = 130 / 800 

Full tree and pruned tree error rates for training data are the same: 0.1625 = 130 / 800

(k)

pred.unpruned <- predict(fit.tree, test, type = "class")
err.unpruned <- round(sum(test$Purchase != pred.unpruned)/length(pred.unpruned), 
    3)

pred.pruned = predict(fit.pruned, test, type = "class")
err.pruned <- round(sum(test$Purchase != pred.pruned)/length(pred.pruned), 3)

Full tree error rate for test data: 0.204.
Pruned tree error rate for test data: 0.204.
Error rate for the pruned tree a little bit increased on the test data.